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Abstract 

We investigate the spatial clustering properties of primordial black holes (PBHs). With minimal 
assumptions, we show that PBHs created in the radiation era are highly clustered. Using the peaks 
theory model of bias, we compute the PBH two-point correlation function and power spectrum. 
For creation from an initially adiabatic power spectrum of perturbations, the PBH power spec- 
trum contains both isocurvature and adiabatic components. The absence of observed isocurvature 
fluctuations today constrains the mass range in which PBHs may serve as dark matter. We briefly 
discuss other consequences of PBH clustering. 
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I. INTRODUCTION 



Primordial black holes (PBHs) are a unique probe of cosmology, general relativity, and 
quantum gravity. Formed by high concentrations of energy density in the early universe, 
PBHs are distinguished from other (astrophysical) black holes by not being formed through 
stellar collapse. In this paper we concentrate on PBHs formed from the direct gravitational 
collapse of density perturbations that are of order unity on the scale of the cosmological 
horizon nn upon horizon entry, though there are other mechanisms for their creation, e. q. 



collapse of cosmic strings |3[ or domain walls [4], or from bubble collisions p in the early 
universe. 

Measurements of the cosmic microwave background (CMB) anisotropy (f| imply that 
density perturbations at the time of decoupling are much smaller (5h ~ 10~ 5 ). As such, PBH 
formation will be cosmologically negligible during and beyond this era. Less constrained are 
the conditions in the early universe before decoupling, and we cannot preclude the existence 
of much larger density contrasts which could have formed PBHs. 

The theory of inflation |2j has been successful in describing both the large-scale homo- 
geneity of the universe and the formation of small-scale structure through the creation of a 
spectrum of cosmological perturbations. It predicts an era of accelerated expansion domi- 
nated by the energy of a slowly rolling scalar field, ending in a period of reheating where 
the energy density is transferred into (more or less) the particles we observe today and the 
radiation dominated epoch begins. The period of reheating is important for PBH produc- 
tion in two ways. First, it is the highest energy scale at which one would expect PBH to 
take place. Gravitational collapse is inhibited by the accelerated expansion, and the number 
density of any PBHs that do form would be drastically diluted. Second, several models of 
inflation exhibit an increase in the amplitude of perturbations at the end of inflation (at the 
epoch of reheating), which increases the probability of PBH formation. 

One topic of interest is the feasibility of PBHs as dark matter (DM) PBHs appear 
to be an a priori good CDM candidate. Formed purely by gravity, they require no special 
extensions to the Standard Model of Particle Physics (such as supersymmetry) , and are 
predicted on quite generic grounds to form in the early universe |2:j. While the smaller 
masses of PBHs (compared to astrophysical black holes) mean that Hawking radiation is 
non- negligible, PBHs that are still in the present universe are still "dark" like other BHs. 
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Because of this, there have been a number of studies of PBHs as CDM in the literature. 
We can split them roughly into three categories: 

QCD PBHs: These are PBHs formed during the QCD phase transition, being a fraction 
of a solar mass 0. Tn is « ^ as evidence from n^cWng even, 

suggested a population of MACHOs in just this mass range. However, in order to produce the 
correct fl m , one needs to invoke a "blue" spectrum {n > 1) of perturbations, which is highly 
disfavored by CMB observations. Further, the evidence that these MACHOS compromise a 
substantial fraction of DM halos is lessening Q| • 

"Spiky" PBHs: These are PBHs formed due to the enhancement of power below a 
certain scale due to features (such as spikes) in the radiation power spectrum 

mm- 

Such a "spiky" power spectrum (a generalizion of a "blue" spectrum, where just the 



power-law slope is changed) can be produced in inflationary models with "plateaus" in the 
inflationary potential. PBHs created in this manner can exist over a larger range of masses, 
given the increased freedom in choosing an inflationary model. Included in this class are 
PBHs created due to perturbation amplification due to preheating 0, 3 [3, Q ■ 

Relic PBHs: These are PBHs of around a Planck mass that exist in some theories 
of quantum gravity as the end result of PBH evaporation j^oi I21I [2^ . As all PBHs with 
initial masses less than ~ 10 15 g would have evaporated by the present day, any model that 
produces a number of light PBHs will leave behind relic PBHs. 

The only limits on PBHs with masses above 10 15 g derive from the requirement that they 
do not overdose the universe (Qpbh < 1), so there is a range of PBH masses over which 
they may serve as DM. Knowing the PBH abundance is necessary, but not sufficient, to fully 
gauge their feasibility as DM. Also important are their spatial clustering properties, as that 
too is constrained by CMB and large scale structure (LSS) data, though to date discussions 
of PBH clustering have been sparse in the literature. A recent general review of PBHs can 
be found in [23j. 



A. PBH Clustering 

The first discussions of PBH clustering came soon after their "discovery" . A theory was 
posited by Meszaros j^J where galaxy formation proceeds from the fluctuations in PBH 
number density. The model does not address how the PBHs are created, but assumes they 
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are around a solar mass and created at or before the QCD phase transition. It claimed that 
for PBH fluctuations that are uncorrected on scales greater than the horizon scale (i.e., 
Poisson fluctuations only), it would be sufficient to able to allow for galaxy formation. This 

n n 

model was refuted in |25( (and later expanded upon in |26(), where it was pointed out that 
the PBH creation process cannot create the "extra" density fluctuations on super-horizon 
scales that was claimed 1 . 

Kotok & Naselsky j3] posit a theory where an initial stage (1st generation) of PBH for- 
mation leads to an early stage of matter (PBH) domination. PBH clustering then enhances 
a second stage (2nd generation) of PBH formation due to collapse in this (pressureless) era; 
specifically, due to the coagulation of PBHs during matter domination. Provided this coagu- 
lation is not complete, the remainder of the 1st generation PBHs evaporate (thus, reheating 
the universe) leaving behind the 2nd generation of PBHs. They claim that with a "blue" 
spectrum of initial perturbations (n > 1.2), PBHs of the 2nd generation are overproduced 
with respect to observational constraints. 

While PBH reheating has been considered it can be shown that that 

the period of PBH domination necessary would lead to the overproduction of (supersym- 
metric) moduli fields and gravitinos upon their evaporation that contradict the predictions 
of big bang nucleosynthesis (BBN). While the authors of 28] seem to confuse the distinction 
between radiation perturbations and PBH perturbations (see their Equation (11)), we show 
later that PBH merging could be a natural consequence of clustering. 

Assuming PBHs comprise the bulk of the CDM, Afshordi, McDonald & Spergel 0] 
study how the discreteness of their population affects the CDM power spectrum. They note 
that PBH perturbations on large scales (super-horizon sized at creation) are a mixture of 
adiabatic (as with other forms of CDM) and isocurvature (due to Poisson fluctuations alone). 
Using Lya forest observations, they use this to constrain the mass of PBHs to be less than 
1O 4 M . They are also the first to investigate PBH cluster dynamics; estimating the lifetime 
due to "evaporation" (different from Hawking evaporation) to show that PBH clusters with 
iV < 3000 objects will evaporate by the current day. We expand on this analysis later. 

Results from microlensing experiments indicate a population of Massive Compact Halo 
Objects (MACHOs) in our galaxy. A possibility that this population is made up of PBHs 



1 Though see [^tJ for a refutation of some of the refutations of 25] and [2^ . 
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of around a half a solar mass, the right mass range for QCD PBHs. In such a population, 



gravitational attraction between PBHs would induce the formation of PBt 



-PBH binaries. 
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38, 



As such, such objects have been studied as sources of gravitational waves 
though to date no such signals have been detected j^J. 

PBHs would be the first gravitationally collapsed objects in the universe. As clustering is 
ubiquitous in other, observed gravitationally collapsed systems (galaxies, clusters of galaxies, 
superclusters, etc), it will be no different for PBHs. The aim of this work is to compute the 
spatial clustering properties of PBHs, and see what impact that has for PBHs in cosmology. 
We will be particularly interested in the viability of PBHs as DM. In Section |H] we describe 
general properties of PBHs we will use throughout the paper. In Section 11111 we derive the 
initial clustering properties of PBHs after their formation, computing the PBH two-point 
correlation function and power spectrum. We conclude in Section IIVI with a discussion of 
observational constraints and avenues for further research. 



II. PBH BASICS 

A black hole of mass M has a Schwarzschild radius Rs = 2GM = j^. Throughout we 
assume that any PBHs have negligible angular momentum and electric charge. 

PBHs form from large perturbations in the radiation density field that are able to over- 
come the resistance of radiation pressure and collapse directly to black holes. For a per- 
turbation of a fixed comoving size, it cannot begin to collapse until it passes within the 
cosmological horizon. The size of a PBH when it forms, therefore, is related to the horizon 
size when the collapsing perturbation enters the horizon 2 . In the radiation dominated regime 
where a oc t 1 ^ 2 and assuming a top-hat window function, the horizon mass is simply 

M H (t) = M P (1) = (2 x 1O 5 M ) (£) (1) 

where tp is the Planck time. Assuming radiation domination, we can rewrite this in terms 
of temperature as 

2 Which is to be expected, being the only characteristic length scale involved. 
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where g* is the effective number of relativistic degrees of freedom. 
The Hubble scale is then determined by the Friedmann equation 



H> = ^P- (3) 



The fluctuation of the (radiation) density field is defined as 



(4) 

and is characterized by its variance on a comoving scale x at a time t as 

a 2 ( X ,t) = ^- 2 J dkk 2 P(k)T(k,t) 2 \W k (k X )\ 2 (5) 

where P(k) is the (primordial) power spectrum, Wk is the Fourier transform of the window 
function, and T(k, t) is the transfer function appropriate for the type of perturbation (adia- 
batic or isocurvature). We further assume that the perturbations are gaussian. It is known 
that the perturbations cannot be completely gaussian, as that would predict perturbations 
with 5 < — 1, implying negative energy densities. This non-gaussianity is especially impor- 
tant in the production of PBHs jlfl| . as they derive from the high end (tail) of the probability 
distribution. Nevertheless, we focus here on the case of underlying gaussian perturbations 
for computational ease. 

Consider a perturbation 8{r H ) smoothed over the comoving Hubble radius r H = R H /a = 
(aH)' 1 . As the underlying perturbations 5 k are assumed gaussian, the smoothed perturba- 
tion will be as well (central limit theorem). As the perturbation must have enough mass to 
overcome pressure, there is a threshold value S c below which a PBH will not form. Further, 
the horizon sized per turbation cannot be larger than unity, or it will pinch off and form a 
separate universe [4l| . Therefore the range that forms PBHs is S G [8 C , 1]. The exact value 
of 5 C is not known precisely. Analytically, 5 C = w , where w is the equation of state parameter 
of the background universe defined through p = wp. The PBH mass was then estimated to 
be M = w 3 / 2 Mh. Numerical simulations of PBH formation 



42, |43j, |44| have shown a more 



complex relation, where 

M = kM h (5 - 5 C )\ (6) 

in accordance with other critical phenomena, where k ~ 3, 7 ~ 0.7 and 5 C ~ 2w. The 
values of these parameters vary depending on the shape of the initial perturbation (gaussian, 
polynomial, etc.). This formally allows for PBHs of an arbitrarily small mass compared 



to the horizon size, though some numerical simulations 42j have showed that there is a 
minimum value of ~ 10~ 3 ' 5 M# as 5 — > 5 C . Rather than focus on one particular formula, we 
can encapsulate our uncertainty in the PBH mass-horizon mass relation with a parameter 
/: 

M PBH = f(w,5 c ,t,...)M H (7) 

While the study of single PBH creation is numerically tractable, the same is not true for 
studying the PBH population as a whole due to their incredible rarity. As such we resort to 
analytical estimates. Given a creation threshold 5 C and the value of the radiation fluctuation 
size at horizon crossing cr ra d(rfr)j the probability of forming a PBH within a given horizon 
volume is simply the probability of having a perturbation with 8 C < 8 < 1, or 

HV-M~ 1/2 ^-^))<« m 

Introducing v = 8 C / <3~rad{ r H) (the threshold in "sigma" units); in the limit where v = 
8 c /a(rH) » 1, the upper limit can be taken to infinity, so that the expression can be 
written in terms of the complementary error function (erfc) as 



a PBB ( V ,M) = e™L = f>(¥-\=JlS=B (.10) 



Note that this can be used to determine the initial PBH density 3 

'PBH _ of M 

pc \M H 

where we have used upbh = P/Vh- 

Without an observed population to compare calculations to, the value of the (physical) 
PBH number density varies in the literature. While j3| did not address PBH formation per 
se, knowing that PBHs form at peaks in the density field implies 

(n + 3) 3 / 2 



n PBH 



(v 2 - l) e~» 2 l 2 R H * (11) 



(2tt) 2 6 3 / 2 

where n is the index of the power spectrum (n = 1 for a scale-invariant spectrum). Whereas 
j3| use peaks in the density field . uses peaks in the metric perturbation to compute a 
density which is identical to the |4a] result, but with (n + 3) replaced with in — 1). This 



3 Sometimes quoted in the literature instead of (3 is a = 13/(1 — 0). In the limit where the PBH mass 
M w Mb, the initial ppbh/ Prad = a. 
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latter calculation only holds for n > 1 47] . Generically, we can write the initial PBH density 

as 

UPBH = (12) 

where N*(y) encapsulates the non-exponential dependence upon v. Equivalently, the initial 
horizon fraction going into PBHs is 

f3 = N,(u)e~ v2/2 . (13) 
The values of NJv) for the different models is summarized below: 



NJv) 



-v 1 erfc approximation 



(14) 



£ (*») ' - 1) BBKS 
l^(¥) 3/ V-D GLMS 
Having determined the initial PBH density, their abundance at subsequent times is simple 
to calculate. PBHs are non-relativistic matter, so Ppbh oc a~ 3 - Because radiation redshifts 
as p r ad oc a" 4 , the PBH to radiation ratio grows until the epoch of matter-radiation equality: 

PPBH,. , B(t) (t*\W 

"a^T * = T^W) ItJ (15) 

After the epoch of equality, Qpbh remains constant during matter domination up until the 
era of vacuum energy domination. The condition that PBHs do not overdose the universe 4 
is Vt PBH (t eq ) < 1/2, or 

1 ( t \ 1/2 

Throughout, we will assume a monochromatic mass function such that Ppbh — Mtipbh- 
We can then write 

B(t) = fN^ 2 ' 2 . (17) 

Figures ^ - El show the lower limit on v derived from the latter equations. This exponential 
dependence of the PBH abundance upon v means we must now then turn to a discussion of 
the form of the underlying power spectrum P(k). 

A given mode crosses within the horizon at a time t given by kn = a(t)H(t). The radiation 
fluctuation on the horizon scale (crossing during radiation domination) is computed using 
Equation (JHJ) 

°Li(r B ,t) = ^J dkk 2 P rad (k)T 2 d (k,t)W 2 (kr H ) (18) 



This is also the condition that PBHs do not induce an early matter-dominated phase. 
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For adiabatic perturbations (which we are assuming for radiation field), T a d(k,t) oc k H 2 (up 
to horizon crossing) and for a power law spectrum P ra d{k) oc k n , 

<ZJTe) « k^ (19) 



48, 



where kn = 1/rjj- The spectrum for n — 1 is known as the Harrison-Zel'dovich 
spectrum (also called a scale-invariant spectrum) and corresponds to fluctuations of different 
physical sizes having identical power when they enter the horizon. Spectra with n > 1 are 
known as "blue" spectra, and correspond to models having more power at smaller scales 
(larger k). 

During radiation domination, the horizon mass Mh oc t oc a 2 oc k]j 2 , or kn oc M^ 1 ^ 2 , so 
that <J 2 ad (r H ) oc M B ~ n ^ 2 . During matter domination the scaling is different, M H ~ pR 3 H oc 
a~ 3 H~ 3 = kjj 3 , or ku oc M H 1//3 , so that a 2 ad oc n ^ 3 . For a pure power law spectrum 
then, we can relate the power at any earlier time to the power today: 

sm- *w> m^m 1 ^ , 2 o, 



vM y \M eq/ 

where subscripts refer to current values and eq refers to the epoch of matter-radiation 
equality. From this, a value of n > 1 can produce sufficient power at small scales to 
produce significant black holes. Our understanding of the physics at these scales in the 
early universe is only theoretical, and thus there may be significant deviations from pure 
power-law behavior then. 



Due to quantum effects 50], a BH of mass M will emit particles as a blackbody with 
temperature T h given by 

T h (M) - — L_ = ML « 10 - (M.Y l e y (21) 
Ih[M) ~ 8nGM~ 8nM i[) \lg) ^ W 

As the temperature is inversely proportional to the mass, this is unobservable for a one solar 
mass (and higher) BH (T/j(M ) pa 62 nK), but cannot be neglected in the mass range of 
PBHs. This emission also corresponds to a mass loss for the PBH, 

•V/ •/-;, -a* SB T^R 2 ) = (22) 

where a* SB is the effective Stefan-Boltzmann constant and is related to the effective number 
of relativistic degrees of freedom in the emitted particles. PBHs therefore have a finite 



9 



lifetime, after which they would have emitted their entire rest mass, given by 

r= ^L_ K(10 -s)(MV (23 ) 

The variation of the parameter a with mass is not great, changing by a factor of 10 over at 
least 7 decades of mass [52]. As the lifetime scales with M 3 , there is a threshold mass above 
which holes will not have evaporated by the present day (to). This threshold mass M* is 
given by 



M* « (4 x 10 14 g 



a(M*) \ / t 



.1/3 



(24) 



Q.9A x 10 25 g 3 /s/ V4.4 x 10 7 s, 
Given the uncertainties in a and t , a threshold mass of M* ~ 10 15 g is typically quoted in 
the literature. 

A large enough abundance of PBHs with M ~ M* will produce a number of observable 
effects through their evaporation in the current day. They would contribute to cosmic rays 
Q.the 7 -ray background 511 keV emission due to positron annihilation in the 
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5611 . Observations 



galactic center |5J] or be the cause of short duration gamma ray bursts 
(or the lack thereof) of PBHs evaporating today depend critically upon not only the number 
density of PBHs present today npBii(to), but also upon how clustered they are within within 
the galaxy. Assuming an isothermal halo model, the effective number density is (npBii{to) 
where ( is the local density enhancement factor Q, H, Q] and ranges from 10 5 — 10 7 . 

PBHs with M < M* would have evaporated by the present day. The main mechanism 
for "observing" PBHs in cosmology is through their Hawking radiation. In the absence of 
a direct detection, the main utility of PBHs is to set limits of PBH abundance at various 
times given a non-detection. Though, PBHs have also been invoked to explain baryogenesis 
H|, reionization and as a solution to the magnetic monopole problem H,H|. 

Evaporating PBHs have their most dramatic effect during the era of BBN, where Hawking 
radiation can alter light element abundances H| . Therefore, the success of BBN implies an 
upper limit to the number of PBHs evaporating at that time. 

Combining Equations (JTJ), ((7j) and (}2l?|) gives the relation 



3a \tp, 

the lifetime r of a PBH created at a time t. What this allows one to do is use information 
from a "late epoch" (time r) to examine conditions at an "early epoch" (time (<t). In 
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'BH abundance from BBN imply 



the above example, r ~ tsBN, and the limits on initial 
P < 1(T 16 for M PBH between 10 9 g and 10 15 g (see, i.e., j7l|). 

This relation depends critically upon the PBH mass monotonically decreasing due to 
evaporation, and not gaining mass in any way (accretion or merging). Should this not be 
the case, the lifetime r is no longer given by the initial PBH mass, and the link between late 
epoch and early epoch is broken. Instead, the energy in PBHs that would have evaporated 



away can now linger for longer periods of time. It was shown in 4l| that PBHs will not 
appreciably increase their mass through radiation accretion. PBH merging then would be 
the dominant mechanism for (significant) mass growth in the radiation dominated epoch. 
Since r oc M 3 , the merging of two equal mass BHs will result in a BH with a lifetime 8 times 
as long. If this merging can continue, then there is a greater chance of PBHs produced in 
the early universe still existing today. 

Depending on the epoch of PBH formation, there is reason to believe there would be 
merging occuring before, say, the epoch of nucleosynthesis, which could skew limits obtained 
from using Equation (f23|) . Assuming an unclustered population, PBH binaries can form in 
the radiation era and be a source of gravitational waves today |22(]. Any PBH clustering will 
only enhance the formation of close PBH binaries (and possibly of larger bound structures), 
and orbital decay will cause merging before evaporation can occur. 



III. BIAS MODEL 



Measuring the two point function (or its Fourier transform, the power spectrum) of as- 
trophysical objects is a powerful tool in studying their clustering properties. The physical 
interpretation of £(r) is as follows. The differential probability of finding two objects (galax- 
ies, clusters, PBHs, etc.) in volume dV\ and (IV2, a distance r apart is given by 

dP = p 2 {l+£(r))dV 1 dV 2 . (26) 

The two point function then measures the excess probability (over random) of finding pairs 
with a separation r (here and throughout we use comoving distances). A large (positive) 
value of £ implies a large amount of clustering (objects are preferentially close to each other), 
a negative value of £ implies anti-clustering (objects are preferentially far away). 

It is important to note that the galaxy-galaxy correlation function £ 99 is not identical 
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to the underlying mass correlation function £ m ; in other words, galaxies are not a perfect 
tracer of mass. Further, different types of objects which may act as tracers (quasar, galaxy 
clusters) have different clustering properties. Measurements of £ for clusters of galaxies 
showed that they were more clustered than galaxies themselves by a factor of 10. Kaiser 
j^l showed that this may be explained using what is now known as the peak-background 
split model of bias: as clusters of galaxies form from higher peaks in the density field than 
galaxies, it is natural that they be more clustered. In the limit of large separation and large 
peaks, the bias is given by 



This can be roughly understood as follows. Split the density field into a long wavelength 
and a short wavelength component. Next, consider a peak in just the long wavelength com- 
ponent ("background"); the physical density field will consist of this component modulated 
by the short wavelength portion. If the threshold for gravitational collapse is close to the 
value of the background peak value, the physical field will cross this threshold a number of 
times in the vicinity of the peak. The regions above threshold, therefore, are preferentially 
found near the background peak. 

The assumptions used are: 

PBH creation is rare: PBH formation occurs during radiation domination (w = 1/3); 
and the radiation perturbations are gaussian. At creation, there will be at most one PBH 
per horizon volume, and PBH formation at around the horizon mass. 

Peaks Theory bias: Since PBH formation is a threshold process, we can use peaks the- 
ory |3] to determine the number density and correlation statistics. While we only consider 
the two-point function and power spectrum here, all higher order correlation functions can 
be derived in a similar manner. 

We now derive the bias for a population of PBHs formed at a single mass scale, compared 
to the underlying radiation field. For the overdensities of PBHs and radiation Spbh and 6 r , 
we define their two point correlation functions 



ipeakir) = — £(r). 



(27) 



£,pbh{t) = (Spbh{x)S P bh{x + r)) 



(28) 



^rad{r) = (S r ad{x)S rad (x + f)) 



(29) 



and the bias parameter 



ipBH{r) = b(r) 2 Crad(r) 



(30) 
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Where, in general, b(r) is not a constant. The averaging done in the definition of the 
correlation functions includes a window function on the scale of the horizon for smoothing. 
Thus, the size of the fluctuations 5 in either the radiation and PBHs is characterized by 

4,o = ^(0) (31) 
From the definition of the radiation and PBH correlation functions, this is given by 

(?PBH,o = b(0)a radfi (32) 
The power spectrum P(k) is defined as 

P(*)=(£)/*r»«r) (33, 

As PBHs form in regions above a certain threshold density, it is straight-forward to 
compute the number density and bias assuming assuming PBHs form at a single mass only. 
The bias is given by an integral over a bivariate gaussian distribution; using the notation of 



Jensen & Szalay 



63j, the full expression is given by 



1 + 




OO POO 



d Vl dy 2 (2n)- 1 [1 - w(r) 



^2^- 1 /2 



(34) 



x exp 



2(1 -w(r) 2 ) 

where w(r) = Crod^V^fado * s the normalized radiation correlation function and v = 
8 c /a ra dfi- It is possible to write this as a power series (the so-called tetrachoric series) 



m w(r\ 



63|, 



oo A 2 



(35) 



where the coefficients are given by 



A, 



yfte^Pertc (^) 



(36) 



5 While the terms perturbation and fluctuation are sometimes used interchangably in the literature to refer to 
an inhomogeneity, we will make a distinction in the usage for radiation and PBHs. The word perturbation 
typically implies smallness in the context of (cosmological) perturbation theory, and we use it to describe 
the (initial) radiation field, as they will be no larger than order unity. As we will show, this will not be 
the case for PBHs, and therefore we use the word fluctuation for their case. 
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where H n are the Hermite polynomials. 



The result of Kaiser 62[ is obtained by assuming w(r) < 1 and v 1, so that only the 
first term in the series need be used to obtain £pRf/(r) ~ v 2 w(r). Relaxing the condition on 
w(r) (but not on u), the coefficients A m —>■ z/ m , obtaining the result of Politzer & Wise ((mJ], 

oo / m\2 



m=l 



-w(r) m = exp(i/ w(r)) 



(37) 



As r 



0, w 



1 by definition, so that in the case of arbitrary u, 



Cpbh(0) 



V^e v2/2 erfc 



^ m-l 
m=l z 



(38) 



Recall that <r% 



'pbho ~ €pbh(0)- F° r large is, it follows that cr 2 DBH0 = e u . In other words, 
PBHs start with a large fluctuation amplitude (compared to radiation) and their evolution 
begins in the nonlinear regime. However, the number density goes as e~ u2 ^ 2 , so the fewer 
PBHs are formed, the more clustered they will be. 

Note this bias is independent of the PBH Mass - Horizon Mass relation (Equation (JTj)). 
Specifically, we have computed the correlation function of horizon sized regions that contain 
at most one PBH. As such PpBH{k) will have an initial upper cutoff at kn- 

While we can compute exactly ^pbh from the peak-background split model, it is custom- 
ary in LSS studies to measure the power spectrum Ppbh instead. Inserting Equation I3TI 
into Equation EH1 we obtain the integral expression 

'sin(Ax)\ 



P 



PBH 



(k) 



4tt 

V 

An r°° 



J drr 2 ^ 



PBH{r 



kr 



(39) 



V 



drr z 



rn 



v 



cxp 



<rad(r) - 1 



radX) 



( sin(Ax) 
y kr 



The lower cutoff at r# = R^/a = k]j , the comoving horizon length at PBH formation, 
is due to the finite size of the PBHs. This will translate into an upper cutoff in Pppuik) 
at kff- Generically, the above integral can be done numerically, but we can say more about 
the nature of the PBH fluctuations without it. 

By expanding the exponential, we can rewrite Equation (J3~9"j) as 



PBH 



(k) 



-Prad(k) + 



radfl 



^ V 

m=2 v 



, o / sin(kr) 
drr 1 



in 



\ kr 







rn 








j m\ 


_ a rad,0 





(40) 
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The higher order terms in the above expansion show the non-linear dependence of Ppbh 
upon P rad . 

Due to the discrete nature of the PBHs, the normalization condition for Ppbh is that as 

k — > 0, Ppbh approaches a spectrum for pure Poisson noise; i.e., a constant value. This is 

manifest in our above expression. The first term, where the PBH power spectrum is simply 

b 2 P ra d, with the bias b given by the Kaiser value of vja. We can Taylor expand the sine 

term in the integrals such that sm(kr)/(kr) — > 1, and those integrals evaluate to constants: 

1 00 4tt u 2m r°° 
Ppo™ = ^E-y-^/ drr%(rr (41) 

V m=2 m[ a rad,0 Jr « 

The total PBH power spectrum then can be written as: 

v 2 

PpBH 

(fc) Ppoisson + ^Prad(k) + Pss(k) (42) 

a rad,0 



U 2 



-irad{r) 



radX) 



(43) 



where 

n /rN 1^~ roc ( (-l) l (kr) 2l \ 4vr 
/=i m=2 h V (2/ + 1)! J ml 
represents the small-scale power when kr is not small. 

To see the behavior of Ppbh at small k, we numerically integrate Equation The 
underlying radiation power spectrum P ra d is a n = 1 spectrum normalized to the four-year 
Cosmic Background Explorer (COBE) value along with a gaussian spike at the horizon scale 
(the latter being normalized to unity). For a fixed 5 C = 2/3, varying the spike amplitude 
will vary the value of v. Figures |U to |H1 show PpBH^k) for four values of v. We see that as 
v increases, the constant (Poisson) power quickly damps out the linear (Kaiser) term. The 
/ = 2 terms of Pss(k) survive for intermediate values of A; as a small negative quadratic 
contribution (oc —k 2 ). 

Note that the power spectrum given in Equation (|4*2*|) is the initial spectrum immediately 
after PBH creation. Due to the different /c-dependence of each of the terms, the power at 
later times (k kn) will not be dominated by the Poisson term. We return to this in 
Section lill CI where we compute the power at horizon crossing at later times. 

From Equation (}3T?j) . we expect Pp i SSO n ~ e^ 2 ; a better fit for v > 4 yields 

Ppoisson ~ y exp Q^ 2 ' 1 ) • (44) 

The power spectrum for a group of iV objects randomly distributed (with a uniform dis- 
tribution) is 1/JV = (nV)' 1 = Note that our above expression for Ppoisson 7^ P , 
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indicating the PBHs are distributed as clusters of objects with mean occupation number 



N c = PpoissonP (45) 

= —N*(v) exp -z/ 2 ' 1 z/ 2 



7 N ' x \4 2 



A. Adiabatic vs. Isocurvature 

We now take an aside to further consider the nature of the PBH fluctuations. That PBHs 
correspond to isocurvature perturbations has been noted in the literature 



35 



45 



though it has not received a lot of attention in the recent PBH publications. In models 
where PBHs constitute the dark matter, it was assumed that their perturbations would be 
purely adiabatic, as with other types of dark matter. We point out that this is not the 
case; a large isocurvature component exists at shorter scales in addition to the adiabatic 
component at longer scales. 

To demonstrate this, assume that radiation is the only component in the universe; there is, 
therefore, no distinction between adiabatic or isocurvature type perturbations. The radiation 
perturbation corresponds to a perturbation in the spatial curvature 6 . Once PBHs are created 
from gravitational collapse, they will evolve as a matter (w = 0) field in the universe. As 
such, we can examine the fluctuations in the PBH density. At the time of PBH creation, 
their fluctuations can be classified as either adiabatic or isocurvature. By assumption, PBHs 
form from the collapse of a density perturbation once it enters the horizon. In the radiation 
dominated era, the PBH mass is close to the horizon mass, so that at most one PBH forms 
per horizon volume. Each PBH is separated by at least a horizon distance. The population 
cannot have correlations on scales smaller than the horizon, so that the perturbations only 
exist for super-horizon scales. Any super-horizon perturbation can be written as a sum of 
adiabatic and isocurvature modes. 

Note that, in our setup, only after PBH creation does the distinction between adiabatic 
and isocurvature perturbations exist. We intend to prove that the PBH fluctuations have 



6 Whether the perturbation is Gaussian or non-Gaussian is largely irrelevent at this point; perturbations 
of order unity must be non-Gaussian to some degree, and we will show in the next section that the 
perturbations of PBHs are generically non-Gaussian. 
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an isocurvature component. This can be generalized to the case where there are additional 
fields and the initial perturbation is wholly adiabatic. 

Isocurvature perturbations correspond to perturbations in the local equation of state 
w — p/p, while adiabatic perturbations correspond to perturbations in the local energy 
density, and thus the local curvature. Consider a volume of space greater than the horizon 
volume at PBH creation. The formation of PBHs cannot change the energy density within 
this space: the gravitational collapse corresponds to a "shuffling" of energy density from 
one form (radiation) into another (matter). The decrement in the radiation energy density 
is exactly balanced by the creation of PBH energy density. Therefore, the curvature is 
unchanged on super-horizon scales. The total perturbation will only become adiabatic if 
this "shuffling" takes place as to satisfy the adiabatic condition. Further, by the second 
law of black hole thermodynamics, a black hole will always have a higher entropy than 
the material that formed it. PBH formation thus corresponds to an increase in entropy, 
and should this process occur non- uniformly, this will result in entropy perturbations, i.e. 
isocurvature perturbations. 

The proof that the PBH fluctuations are isocurvature, then, derives from the fact that 
PBH formation is highly non-uniform. Equivalently, that PBHs are created highly clustered, 
which was shown earlier in this section. Using the notation of |7j], we write the entropy 
perturbation as 



where 5' rad is the radiation perturbation after PBH formation, and S' rad ^ S ra d. Using the 
parameter B from Equation (jl(J|) . we can trivially write 



SpBH — 5 PBH — ~;Kad 



(46) 




(47) 



which allows us to relate the perturbations as 




(48) 



We can use this latter equation to rewrite Equation (|47)|) as 




(49) 



&PBH — -T^rad 



(50) 
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which is now a function of the initial radiation perturbation and the (final) PBH fluctuation, 
and the approximation holds as long as B <C 1. While 5 ra d < 1 by assumption, we know 
from Equation that Spbh typically will not due to clustering. We see that the entropy 
perturbation is simply a function of the bias parameter: 



It is apparent that the isocurvature perturbation is almost inevitable for realistic (rare) PBH 
production. The bias parameter b will be dependent on scale; in Fourier space b is given 



roughly by J Ppbh / Prad- For a given k, the bias is dominated by the term domination the 
power spectrum as given in Equation At the smallest scales (close to PBH creation 

scales), the bias is largest and using Equation (|37j) gives b ~ exp(i/ 2 /2)/a 3> 1. For larger 
scales, the linear (Kaiser) bias gives b = vjo. In either case, the parameters (a, v) would 
have to be finely tuned in order to produce a purely adiabatic PBH perturbation. 

We note that this mechanism for generating an isocurvature perturbation is independent 
of the process that created the initial (adiabatic) perturbation, though we assume through- 
out that it is done through an epoch of cosmological inflation. This mechanism then is 
an exception to the generally held thought that an isocurvature perturbation cannot be 
produced from single field inflation The reason this occurs is that PBH creation {i.e. 
gravitational collapse) is an inherently non-linear and non-perturbative process that is not 
bound by this restriction from perturbation theory. PBH dark matter is not like particulate 
dark matter. Further, for PBHs lighter than M* this isocurvature fluctuation is transferred 
to the products of Hawking evaporation. Thus, the absence of an observed isocurvature 
perturbation implies a limit on the number of PBHs that have evaporated in the past. We 
plan to further explore this topic in a future paper. 

B. Gaussian vs. Non-Gaussian 

In our derivation of the PBH number density and clustering properties, we assumed the 
underlying radiation perturbation was gaussian. As PBHs form only at the peaks of the 
density field, and the initial size of the fluctuation is greater than unity, the PBH fluctuations 
cannot be gaussian. They appear instead to be lognormal (LN) in character. Roughly, a 
LN distribution is the exponentiation of a gaussian distribution. The two-point correlation 




(51) 
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function of a LN field is given by |79( 

l + ^(r) = exp(S(r)), (52) 

where H(r) is the correlation function of a gaussian field with variance H(0) = S 2 . From 
Equation (J37j) . this is exactly the correlation function for the PBH population assuming 
S(r) = v 2 w(r) and ^ ^> 1. 

An isocurvature perturbation necessarily is defined between two components, here radi- 
ation and PBHs. While we have been focusing on the PBHs, there is of course a change 
in the radiation perturbations; the increase in PBH density is exactly cancelled by a de- 
crease in radiation density. From the perspective of the radiation field, not only is there 
an isocurvature component along with the (initial) adiabatic component, but there is now 
a non-gaussian fluctuation along with the (initial) gaussian one. 



C. PBH Fluctuation evolution 



The evolution of the PBH population after creation is a complex problem, outside the 
bounds of perturbation theory due to the size of the initial PBH fluctuations, and better 
addressed as an N-body problem However, we can make a rough estimate of the power 
at horizon crossing of other scales using the results from cosmological perturbation theory. 
We can break the PBH power spectrum into its isocurvature and adiabatic components: 

PpBH(k) = (Pp BH {k) - ^Pradik)^ + ^P rad (k) 

= P iso {k) + P ad {k). (53) 

We can then write the variance at horizon crossing as 

v 2 PB H (r H ,t) = ^- 2 J dkk 2 (P lso (k)T 2 so (k,t) + P ad (k)T 2 ad (k,t)) W 2 (kr H ). (54) 

Rather than computing this explicitly, we will note that for power law spectra where 
Pi SO (k) oc k Uiso and P a d{k) oc k n , their contributions to the variance at horizon crossing can 
be written as 

1— n 1—n 

•K-> - -iwrt (£) 1 . m 

2 ,TT-U ( M e q \~^ M H \ 2 



\M J \M, 



(56) 
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while the total variance is their sum: 



a 



2 ~ °lo + °ld- (57) 



The condition for scale-invariance is no scaling with mass; for adiabatic perturbations 
this is n = 1, for isocurvature perturbations this is rii SO = —3. While the adiabatic portion 
of PpBii{k) is scale- invariant by assumption, for scales larger than the horizon size at their 
creation, the isocurvature component has a flat spectrum (n, so ~ 0) and diminishes at longer 
scales. Thus while the isocurvature portion dominates initially, there is a crossover scale 
where the spectrum becomes adiabatic. Given the lack of measured isocurvature component 
at the time of the CMB (upper limit on isocurvature fraction is fi SO < 0.33, from 8J|), we 
can put a limit on the PBH population so that it does not violate this bound. Roughly, at 
the scale of matter-radiation equality (M E q ~ 10 48 g), 



° (r EQ ) =af so + a 2 ad = 5 2 H} (58) 



and the bound is 



^isoiTEQ) < fLSl, ( 59 ) 

where 5h = 1.91 x 10 -5 . To compute crf so (r E Q), we assume Pi SO {k) ~ Ppoisson, which, as 
shown in Figures |U - [7[ is valid for k < fc^/lO. The upper limit is plotted in Figures [T] - El 
for three different values of /. This constraint becomes an upper limits on the (initial) PBH 
mass if it is to serve as dark matter. For / = 1, allowed regions for PBH dark matter all 
have Mpbh < M Q , so that there would be no confusion with astrophysical BHs. As we 
decrease /, the upper limit increases: for / = 10~ 3 ' 5 it is pushed above the confusion limit. 



IV. CONCLUSIONS 



We have shown that for PBHs to serve as dark matter, clustering constrains them to lie 
in a particular mass range. Further, PBHs will preferentially be found in clusters. 

As shown in the previous section, PBH fluctuations enter the horizon with a very large 
amplitude {opbh ~ e^ 2 ). It is therefore no longer value to treat their evolution using linear 
perturbation theory, as one is able to do for other forms of CDM. Instead, we examine the 
sub-horizon evolution of the PBH population as an N-body problem. Being non-relativistic, 
PBHs will cluster hierarchically (just as CDM); creating smaller bound systems that get 
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incorporated into larger ones. The internal dynamics of these systems are determined solely 
by gravitational clustering, analogous to other gravitationally bound systems such as star 
clusters and galaxies. For this, we are aided by the work done in the context of studying more 



massive black holes in globular clusters 



65j | and galaxies |66j]. In those cases, gravitational 



interactions tend to either produce bound pairs or ejections, rather than BH coalescence 



What occurs in the case of PBHs depends upon how many form in a "PBH cluster" and 
what their initial separations are. The estimate of cluster population in Equation is 
likely an overestimate since our approximation for S,pbh{^) breaks down for small r. The 
initial separations should be on order the horizon size at formation, being the only length 
scale involved in PBH formation. This would seem to indicate rather compact clusters (initial 
separation on order the size of the PBHs themselves), though more work (e.g., higher order 
statistics, numerical simulations) is needed to verify this. 

Frequent merging due to clustering could have a profound impact upon cosmology. Since 
their lifetime r oc M 3 , PBHs, due to merging, exist longer than they would have initially. 
This could feasibly lead to a PBH population in the present universe that was formed in the 
earliest moments of the early universe, opening up a new and unique observational window 
into that time. At the very least, PBH merging in clusters dramatically changes the limits 
on initial PBH abundance, such as those used to put limits on models of inflation 



74, 



75 



69 



^fij . The issue of PBH clusters and merging is discussed more fully in 

a companion paper 

Limits on the current number density of PBHs depend critically upon how clustered PBHs 
are in our galaxy. Naively, from our work in this paper, we might expect a local clustering 
enhancement £ ~ e 1 ' 2 / 2 , or £ ~ 10 22 for v = 10. This is many orders of magnitude larger 
than the factors of 10 7 computed in the literature. This ignores the effect of PBH merging 
though; sufficient merging might concentrate all galactic PBHs into the center SMBH. This 
will have implications for models where PBHs are used to be the "seed" BHs needed for the 
growth of SMBHs in the centers of galaxies 

This PBH merging scenario we have discussed has other predictions. One prediction is 
more gravitational wave emission than originally assumed for a uniform PBH population. 
This is due to the increased probability of PBH binary formation and emission from resonant 
bound states. 
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The other prediction is related to dark matter. Suppose now that PBHs are not the 
only component of the dark matter, and that there also exists a "standard" CDM candidate 
with adiabatic perturbations (in accordance with CMB measurements), in which case the 
CDM perturbation amplitude is related to the radiation perturbation amplitude by Scdm = 
(3/4)<W 

Perturbations in the radiation density can only collapse (into PBHs) if they are of suffi- 
cient amplitude on the scale of the horizon. Perturbations smaller than this, in accordance 
with linear perturbation theory, will simply oscillate, but not collapse. This implies that 
there will be scales slightly larger than those where PBH formation took place where S r is 
below the threshold for PBH formation but still large compared to, say, the amplitude at 
the time of the CMB (10~ 5 ). There is, accordingly, a similarly large perturbation in the 
CDM density assuming adiabaticity. While the linear growth of matter perturbations is 
delayed until after matter-radiation equality, they still grow logarithmically in the radiation 
dominated era. This leads to the possibility that they will become non-linear before equality, 
and forming bound dark matter structures along with PBHs. In which case, one would have 
to include the interaction between these two populations of primordial bound objects. 
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FIG. 1: Allowed region in v — MpBH parameter space for PBHs, assuming / = 1. Solid curve is the 
upper limit on v due to isocurvature perturbations (from Equation 1)59(1 ). Other curves are lower 
limits on v due to number density (Equation ((14(0 : long dashed line uses the erfc approximation, 
dotted line uses the BBKS formula with n = 1, short dashed line uses the GLMS formula with 
n = 1.5. Heavy lines show where PBH dark matter is allowed by the isocurvature constraint. 
Shown also is the temperature of the universe T when PBHs form. The line at M = M* ~ 10 15 g 
is mass below which PBHs would have Hawking evaporated by the current day (assuming no 
accretion or merging). The line at M ~ 3Mq is the mass above which PBHs would be confused 
with astrophysical BHs. 
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FIG. 2: The same as Figure^ except with / = 0.1. 
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FIG. 4: The PBH Power Spectrum for v = 1.17. Dotted line is the radiation power spectrum, 
consisting of a n = 1 spectrum with COBE normalization, along with a gaussian spike in power 
at k = 1. Solid line is the PBH power spectrum, dashed line is the quadratic estimate of the PBH 
power spectrum. 
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FIG. 5: The PBH Power Spectrum for v = 2.62. Dotted line is the radiation power spectrum, 
consisting of a n = 1 spectrum with COBE normalization, along with a gaussian spike in power 
at k = 1. Solid line is the PBH power spectrum, dashed line is the quadratic estimate of the PBH 
power spectrum. 
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FIG. 6: The PBH Power Spectrum for v = 3.71. Dotted line is the radiation power spectrum, 
consisting of a n = 1 spectrum with COBE normalization, along with a gaussian spike in power 
at k = 1. Solid line is the PBH power spectrum, dashed line is the quadratic estimate of the PBH 
power spectrum. 
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FIG. 7: The PBH Power Spectrum for v = 8.30. Dotted line is the radiation power spectrum, 
consisting of a n = 1 spectrum with COBE normalization, along with a gaussian spike in power 
at k = 1. Solid line is the PBH power spectrum, dashed line is the quadratic estimate of the PBH 
power spectrum. 
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